function loadaggdata(SampleStart,SampleEnd)

    agg_data = readdlm(dataDir * "aggvar_sim.csv", ',', Float64, '\n'; skipstart=0);
    return agg_data[SampleStart:SampleEnd,:]

end

##

function loadaltdata(SampleStart,SampleEnd)

    # dataDir  = "$(pwd())/data/Para" * nDataSel *"/"
    percentiles_data     = readdlm(dataDir * "percentiles_data_fVAR3.csv", ',', Float64, '\n'; skipstart=1)
    gini_coef_data       = readdlm(dataDir * "gini_coef_data_fVAR3.csv", ',', Float64, '\n'; skipstart=1)
    probMass_aboveX_data = readdlm(dataDir * "probMass_aboveX_data_fVAR3.csv", ',', Float64, '\n'; skipstart=1)

    percentiles_data     = percentiles_data[SampleStart:SampleEnd,2:end] # include 10, 20, 50, 80, 90th percentiles
    vec_percs            = [0.1; 0.2; 0.5; 0.8; 0.9]
    vec_cutoffs          = [1.0; 1.5; 2.0]
    probMass_above15_data = probMass_aboveX_data[SampleStart:SampleEnd,3]
    probMass_above20_data = probMass_aboveX_data[SampleStart:SampleEnd,4]
    gini_coef_data       = gini_coef_data[SampleStart:SampleEnd,2] # include gini coefficient

    return percentiles_data, vec_percs, gini_coef_data, probMass_above15_data, probMass_above20_data

end

##

function loaddensdata(SampleStart,SampleEnd,K,nfVARSpec)

    sNameFile = "K" * string(K) * "_fVAR" * nfVARSpec
    PhatDensCoef_factor = readdlm(loaddir * sNameFile * "_PhatDensCoef_factor.csv", ',', Float64, '\n'; skipstart = 1);
    PhatDensCoef_lambda = readdlm(loaddir * sNameFile * "_PhatDensCoef_lambda.csv", ',', Float64, '\n'; skipstart = 1);
    PhatDensCoef_mean   = readdlm(loaddir * sNameFile * "_PhatDensCoef_mean.csv", ',', Float64, '\n'; skipstart = 1);
    MDD_GoF             = readdlm(loaddir * sNameFile * "_MDD_GoF.csv", ',', Float64, '\n'; skipstart = 1);
    Vinv_all            = readdlm(loaddir * sNameFile * "_Vinv_all.csv", ',', Float64, '\n'; skipstart = 1);
    N_all               = readdlm(loaddir * sNameFile * "_N_all.csv", ',', Float64, '\n'; skipstart = 1);

    PhatDensCoef_factor = PhatDensCoef_factor[SampleStart:SampleEnd,:]
    MDD_GoF             = MDD_GoF[SampleStart:SampleEnd]
    Ktilde              = size(PhatDensCoef_lambda)[1]
    SampleT             = SampleEnd-SampleStart+1

    # Load ME covariance matrices
    # Needs to be adjusted for Lambda
    VinvLam_all = zeros(Ktilde, Ktilde, SampleT)
    tt_sel = 1
    for tt = 1:size(N_all)[1]
        if SampleStart <= tt <= SampleEnd
            Vinv_t = Symmetric(Vinv_all[K*(tt-1)+1:K*tt,:])
            VinvLam_all[:,:,tt_sel] = PhatDensCoef_lambda*Vinv_t*PhatDensCoef_lambda'*N_all[tt]
            tt_sel = tt_sel + 1
        end
    end

    return PhatDensCoef_factor, MDD_GoF, VinvLam_all, PhatDensCoef_lambda, PhatDensCoef_mean

end

# function loadknots(quant_vec, quant_sel)
#
#     dataDir = "$(pwd())/data/"
#     densdraws_data   = CSV.read(dataDir * "simul_asset_data_final_Tlong.csv", DataFrame, header = false);
#     employdraws_data = CSV.read(dataDir * "simul_employ_data_final_Tlong.csv", DataFrame, header = false);
#
#     #simul_data  = convert(Array,simul_data)
#     densdraws   = convert(Array,densdraws_data)
#     employdraws = convert(Array,employdraws_data)
#
#     # subsequently use the same knots regardless of sample size N,T
#     knots_all = quantile(densdraws[employdraws.==1], quant_vec)
#
#     ii=getindex.(findall(K_vec.-K.==0),[1 2])[1] # find index ii where K==K_vec
#     knots = knots_all[quant_sel[ii,:].==1]
#
#    return knots
#
# end
